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ABSTRACT  .^Procedures  are  developed  for  constructing  asymptotic  solutions 
for  certain  nonlinear  singularly-perturbed  vector  two-point  boundary  value 
problems  having  boundary  layers  at  one  or  boLh  end  points.  The  asymptotic 
approximations  are  generated  numerically  and  can  either  be  used  as  is  or  to 
furnish  a  two-point  boundary  value  code  (e.g.  COLSYS)  with  an  initial  approx¬ 
imation  and  a  nonuniform  computational  mesh.  The  procedures  are  applied  to 
several  examples  involving  the  deformation  of  nonlinear  elastic  beams. ^ _ 

1.  INTRODUCTION.  We  consider  singularly-perturbed  two-point  boundat-y^ 
value  problems  for  nonlinear  vector  systems  of  the  form 

x  -  f ( x , y , t , c)  ,  cv  -  g(x,y,t,e)  ,  0  <  t  <  1  (la,b) 

a(x(0),y(0),e)  =  0  ,  b(x( 1) ,y( 1) , c)  -  0  (lc.d) 

where  x,  y,  a,  and  b  are  vectors  of  dimension  m,  n,  q,  and  r  *  m  +  n  -  q, 
respectively.  We  seek  to  find  limiting  solutions  of  problem  (1)  as  the  small 
positive  parameter  e  tends  to  zero;  however,  to  do  this  in  complete  generality 
is  very  difficult  and  beyond  the  grasp  of  our  current  understanding.  Ulus,  we 
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simplify  problem  (1)  considerably  by  assuming,  in  addition  to  natural 
smoothness  hypothesis,  that  (i)  g,  a,  and  b  are  linear  functions  of  the  fast 
variable  y,  i.e. 

g(x,y,t,e)  -  g](x,t,e)  +  G2(x,t,e)y  (2a) 

a(x(0) ,y(0) ,e)  -  aj(x(0) , t)  +  A2(x(0) , c)y(O)  (2b) 

b(x( 1) ,y(l) ,e)  =  bj(x( 1) , e)  +  B2(x( 1) , e)y( 1)  (2c) 

(ii)  that  G2(x,t,e)  has  a  hyperbolic  splitting  with  k  >  0  stable  eigenvalues 
and  n  -  k  >  0  unstable  eigenvalues  for  all  x  and  0  <  t  <  1,  and  (iii)  that  q  > 
k  and  r  >  n  -  k. 

With  the  assumed  hyperbolic  splitting,  we  would  expect  y  to  vary  rapidly 
relative  to  (the  slow  vector)  x  in  narrow  boundary  layer  regions  near  boLh  t  = 
0  and  1.  We  thus  seek  limiting  solutions  having  the  form 

x( t, e)  =  X(t)  +  0(e)  ,  y(t,E)  =  Y(t)  +  p(  t)  +  v(o)  +  0(e)  (3a, b) 

where  the  initial  layer  correction  u(t)  and  the  terminal  layer  correction 
v(o) ,  respectively,  decay  to  zero  as  the  stretched  variable 

t  =  t/c  or  u  -  ( 1  —  t )  / 1.  -  (4a, b) 

tend  to  infinity.  The  limiting  solution  X(t),  Y(t)  within  0  <  t  <  1  must 
necessarily  satisfy  the  reduced  system 

X  «=  f (X,Y,t ,0)  ,  0  -  g( X,Y ,t ,0)  (3a, b) 

Because  G2  is  everywhere  nonsingular,  we  can  use  Eqs.  (2a)  and  (5b)  to 
determine 

Y(t)  =  -G2_I(X,t,0)gi(X,t,0)  (6) 

in  a  locally  unique  way,  and  there  remains  the  mth  order  nonlinear  differen¬ 
tial  system  (Eq.  (5a))  for  determining  X(t). 

In  order  to  completely  specify  the  reduced  solution  we  must  prescribe  m 
boundary  conditions  for  equations  (5a) .  We  do  this  by  providing  a 
"cancellation  law"  which  selects  a  combination  of  q-k  initial  conditions  (Eq. 
(2b))  and  of  r  -  n  +  k  terminal  conditions  (Eq.  (2c))  to  be  satisfied  by  X  and 
Y.  In  Section  2  we  present  a  numerical  procedure  for  determining  the  boundary 
conditions  for  the  reduced  system  that  uses  an  orthogonal  iratrix  E(x,t)  to 
reduce  the  matrix  G2(X(t),t,0)  to  a  block  Iridiagonal  form  so  that  the  stable 
and  unstable  eigenspaces  may  be  separated.  The  boundary  layer  corrections 
u(t)  and  v(o)  in  Eqs.  (3)  compensate  for  the  cancelled  initial  and  terminal 
conditions,  respectively,  and  they  can  be  determined  once  X(t)  has  beer. 


In  Section  3  we  discuss  a  numerical  procedure  for  determining  the 
asymptotic  approximation  (Eq.  (3))  which  uses  the  general  purpose  two-point 
boundary  value  code  COLSYS  to  solve  the  reduced  problem  and  then  adds  numeri¬ 
cal  approximations  to  the  boundary  layer  corrections.  This  approximation  is 
considerably  less  expensive  to  obtain  than  solving  the  full  stiff  problem 
numerically  and  it  has  the  advantage  of  improving  in  accuracy,  without  any 
additional  computational  cost,  as  the  small  parameter  c  tends  to  zero.  How¬ 
ever,  when  e  is  only  moderately  small  our  asymptotic  approximation  may  not  be 
sufficiently  accurate  for  some  purposes,  so  we  have  developed  a  procedure  (cf. 
Section  3)  that  generates  an  improved  solution  by  using  COLSYS  to  solve  the 
complete  problem  (Eqs.  (1)  and  (2))  with  our  asymptotic  approximation  as  an 
initial  guess.  In  order  for  this  approach  to  succeed  we  must  also  provide 
COLSYS  with  an  initial  nor.uniform  mesh  Lhat  is  appropriately  graded  in  Lhe 
boundary  layers  (cf.  Ascher  and  Weiss  (Per.  2))  and  we  give  an  algorithm  for 
constructing  such  a  mesh  in  Section  3.  While  ;>ur  procedure  d"es  not  appear  to 
be  optimal,  we  show  by  an  example  involving  the  deformation  cf  a  nonlinear 
elastic  beam  (cf.  Section  4)  that  it  does  offer  some  advantage  over  Lhe  more 
standard  approach  of  continuation  in  c,  where  one  starts  with  a  large  value  of 
e  (e.g.  e  »  1)  and  a  crude  initial  guess  and  reduces  e  in  steps  so  that  the 
mesh  is  gradually  concentrated  into  boundary  layer  regions. 

We  close  Section  4  with  a  second  nonlinear  beam  example  that  is  beyond 
Lhe  capabilities  of  our  present  methods  because  the  matrix  C2  is  a  function  of 
y.  Flaherty  and  O'Malley  (Ref.  6)  analyzed  this  problem  and  showed  that  its 
solution  becomes  unbounded  as  £  *  0.  We  include  the  numerical  solution  of 
this  problem  in  this  paper  in  order  to  show  one  of  the  mar.v  challenging 
effects  that  can  occur  with  singularly-perturbed  problems. 

Finally,  in  Section  5  we  discuss  cur  results  and  present  some  suggestions 
for  future  investigations. 

2.  ASYMPTOTIC  APPROXIMATION.  In  order  to  calculate  Lhe  houndary 
conditions  for  the  reduced  problem  (Eqs.  (5a)  and  (6))  and  Lhe  boundary  layer 
corrections  u(t)  and  v(o)  we  calculate  Lhe  Schur  decomposition  of  the  matrix 
G2  at  t  -  0  and  L  =  1.  In  particular,  at  t  =  0  we  find  ar.  orthogonal  matrix 
E(x(0))  such  that 

T_(x(0) )  U(x(Q) ) 

G2(x(0) ,0,0)E(x(0))  «  E(x( 0) )  (7) 

0  T+(x(0) ) 


where  T_  is  k  x  k  and  upper  triangular  with  Lhe  stable  eigenvalues  of  G2 ,  and 
T+  is  upper  triangular  with  Lhe  r.-k  unstable  eigenvalues  of  G2.  The 
decomposition  (Eq.  (7))  can  often  be  obtained  analytically;  however,  when  this 
is  not  possible  or  practical  it  can  be  determined  numerically  by  using  the  QR 
algorithm  (cf.  Golub  and  Wilkinson  (Ref.  7)  and  Ruhe  (Ref.  9)  for  specific 
procedures) . 


We  partition  E  after  its  kth  column  as 

E  =*  [  E_  E_]  (8) 

and  note  that  E_  spans  the  stable  eigenspace  of  G2  at  t  *  0  and 

P  -  E_  E_t  (9) 

is  a  projection  onto  this  eigenspace. 

Near  t  -  0,  we  assume  that  the  terminal  layer  correction  v  is  negligible, 
substitute  the  asymptotic  approximation  (Eq.  (3))  into  the  differential 
equations  vEqs.  (la,b)),  use  the  reduced  system  (Eq.  (5)),  and  retain  only  the 
leading  order  terms  to  find  that  u(t)  satisfies  the  conditionally  stable 
system 


dy/d t  -  G2(0)y  (10) 

where  (here  and  below)  we  use  the  argument  t  to  denote  conditions  evaluated  at 
x(t)  *  X(t),  t,  and  e  ■»  0,  e.g., 

G2(0)  :=  G2(X(0) ,0,0)  •  (11) 

Integrating  Eq.  (10) 

G2(0)t 

p(t)  -  e  y(0)  (12) 

We  require  that  y(r)  decays  as  t  Increases  and  this  will  be  the  case  provided 
that  y(0)  is  in  the  stable  eigenspace  of  G2(0);  thus,  using  Eq.  (9)  we  require 

U(0)  -  P(0)m(0)  -  E_(0)E_t(0)m(0)  (13) 

Using  Eqs .  (3),  (13),  and  (2b)  in  Eq.  (lb)  we  find  that  the  limiting 
initial  conditions  have  the  form 

ai(0)  +  A2(0)  (Y(0)  +  L-(0)E_t(0)w(0)]  -  0  (14) 

We  assume  that  A2(0)E_(0)  has  its  maximal  rank  k  and  construct  a  q  x  q  matrix 

lT  =■  [L-T  U-T]  (15a) 

that  reduces  it  to  row  echelon  form,  i.e., 


—  -  — 

L_ 

V_ 

- 

A2(0)E_(0)  = 

L_ 

0 

_ 

(15b) 


where  V_  is  k  x  k  and  nonsingular.  Multiplying  Eq.  (14)  by  L  and  using  Eqs. 
(13)  and  (15)  gives  the  initial  layer  jump  and  the  q-k  initial  conditions  for 
the  reduced  problem,  respectively,  as 


11(0)  -  -E_(0)V_-1L.[ai(X(0)  ,0)  +  A2(X(0)  ,0)Y(0) ) 


(16a) 


and 

f(X(0))  :«  1— ( aj( X(0) ,0)  +  A2(X(0) ,0)Y(0) )  -  0  .  (16b) 


We  find  the  terminal  layer  jump  and  the  r  -  (r.-k)  terminal  conditions  for 
the  reduced  problem  in  an  analogous  fashion  with  the  exception  that  we  define 
E(x(l))  such  that 


G2(x( 1) ,1,0) E(x( 1) )  -  E(x( 1 ) ) 


T+(x( 1) ) 
0 


U(x(  1)) 
T-(x(l)) 


(17) 


which  we  partition  after  its  (n-k)th  column 

E  -  [E+  E+] 


(18) 


In  parallel  with  Eqs.  (7)  and  (8),  the  matrices  T_,  T+,  and  E+  contain  the  k 
stable  eigenvalues,  the  n-k  unstable  eigenvalues,  and  span  the  unstable  eigen- 
space,  respectively,  of  G2  at  t  *  1.  Our  reasons  for  switching  the  positions 
of  the  matrices  containing  Lhe  stable  and  unstable  eigenvalues  of  G2  Is  that 
there  is  no  simple  and  stable  computational  procedure  for  finding  a  set  of 
vectors  that  span  a  given  subspace  and  are  noi  in  Lhe  leading  columns  of  an 
orthogonal  matrix  like  E  (cf.  Golub  and  Wilkinson  (Ref.  7)). 


Now,  following  the  procedure  that  we  used  for  the  initial  layer,  we  find 
that  the  terminal  layer  correction  satisfies 


G2(1)o 

v(o)  »  e  v(  0)  (19) 

In  order  for  v(o)  to  decay  as  o  increases  we  require  v(o)  to  be  in  the 
unstable  eiger.space  of  G2(l);  thus,  we  take 

v(0)  -  QUMO)  -  M  1)E+T(1)  v(0)  (20) 

where  Q  is  a  projection  onto  the  (n— k)  dimensional  unstable  eieenspace  of 
G2(l). 

We  assume  that  B2(1)E+.(1)  has  its  maximal  rank  n-k  and  find  a  r  x  r 
matrix 


RT  -  [k+t  k+t] 

that  reduces  it  to  the  row  echelon  form 


(21a) 


b2(  DM  l)  - 

v+ 

R+ 

0 

(21b) 


1 


1 

where  V+  is  (n-k)  x  (n-k)  and  nonsingular.  Multiplying  Eq.  (Id)  by  R,  using  i 

Eqa.  (2c),  (3),  (20),  and  (21),  and  retaining  only  the  leading  order  terms  we 
find  the  terminal  layer  jump  and  the  r  -(n-k)  terminal  conditions  for  the 
reduced  problem,  respectively,  as 

v(0)  -  -E+tDV+^R+tb^XO)  ,0)  +  B2(X(  1)  ,0)Y(  1)  )  (2*2a) 

and 

*(X(1))  R^lb^Xd)  ,0)  +  B2(X(1),0)Y(1)]  -  0  (22b) 

In  the  interest  of  brevity,  we  have  omitted  several  details  of  our 
construction  and  have  not  attempted  to  justifv  the  asymptotic  validity  of  our 
procedure.  These  topics  will  be  the  subject  of  a  forthcomi n.»  paper  by 
O'Malley  and  Flaherty  (kef.  8). 

3.  NUMERICAL  PROCEDURE.  Our  computational  procedure  consists  of  first 
solving  the  reduced  problem  (cf.  Eqs.  (5a),  (6),  (16b),  and  (22b))  numerically 
and  then  adding  any  boundary  layer  corrections.  Since  the  reduced  problem  is 
not  stiff  we  can  use  any  good  code  for  two-poinL  boundary  value  problems  (cf. 

Childs  et  al.  (Ref.  3))  and  we  have  chosen  to  use  the  collocation  code  COLSYS 
of  Ascher,  Christiansen,  and  Russell  (Ref.  1). 

Since  the  reduced  problem  is  generally  nonlinear  and-  since  COLSYS  solves 
nonlinear  problems  using  a  damped  Newton  method  we  have  Lo  supply  formulas  for 
evaluating  the  Jacobians  of  f,  Y,  $,  and  Y  with  respect  to  X.  We  do  this  by 
providing  analytical  formulas  for  these  Jacobians  that  neglect  the  influence 
of  the  derivatives  of  E,  L,  R,  and  G2.  This  procedure  has  not  failed  on  any 
of  our  examples;  however,  an  alternate  possibility  would  be  to  approximate  the 
Jacobians  by  finite  differences. 

We  start  the  Newton  iteration  with  a  uniform  mesh  and  the  default  initial 
guess  X<0)(t)  for  X(t)  that  is  provided  by  COLSYS  and  calculate  successive 
approximations  X^P^(t)  until  convergence  is  attained.  AL  each  iteration  step 
we  calculate  an  approximation  E^P^(t)  to  E(t)  for  t  =  0  and  1  as  the  Schur 
decomposition  of  G2(X^PHt)  ,t,0)  .  In  the  examples  of  Section  4  we  used 
analytical  formulas  for  E  rather  than  the  numerical  procedures  of  Golub  and 
Wilkinson  (Ref.  7)  or  Ruhe  (Ref.  9).  Finally.  L^P^  and  r(p)  are  obtained 
using  Gaussian  elimination  to  row  reduce  A2( X'  P)(0)  ,0)  E_(p)(0)  and 
B2(X?P>(1)  ,0)E*(P>(1) ,  respectively. 

When  the  above  procedure  converges  we  calculate  boundary  layer 
corrections  p(t)  and  v(o),  for  a  given  value  of  e,  using  Eqs.  02),  (16a), 

(19),  and  (22a),  and  add  these  to  the  reduced  solution  in  order  to  get  the 
0(e)  asymptotic  approximation  (Eq.  (3)).  For  moderately  small  values  of  e 
this  approximation  may  not  provide  a  sufficiently  accurate  representation  of 
the  solution  and,  in  this  case,  we  use  it  as  an  initial  guess  to  COLSYS  and 
solve  the  complete  problem  (Eq.  (1)).  Unfortunately,  this  procedure  will  fail 
unless  we  also  provide  COLSYS  with  an  initial  nonuni  form  partition 


it 


{0  -  t0  <  ti  < 


•  •  • 


<  tN  -  1} 


(23) 


that  is  appropriately  graded  within  the  boundary  layers.  We  seek  to  find  u  so 
that  the  poir.lwise  error  satisfies 


|  I e C t-i )  |  |  <  6(1  +  |  | u( 1 1 )  |  | )  ,  i  -  1,2,...,N-1  (24) 

where  6  is  a  prescribed  tolerance,  uT  [xT,yT] t  e  is  the  difference  betwfeen 
u  and  its  collocation  approximation,  and 

| lu(ti) | )  max  | u j ( t i ) j  (25) 

1  < j  <ra+n 

We  have  based  our  condition  for  determining  it  on  a  pointwise  error  criteria 
since  this  seemed  to  work  better  in  practice  than  a  global  criteria.  This  is 
somewhat  surprising  since  COLSYS  uses  a  global  error  criteria  to  select  a 
mesh. 


We  assume  that  the  final  partition  selected  by  COLSYS  to  solve  the 
reduced  problem  satisfies  equation  (24)  outside  of  boundary  layer  regions  and 
we  seek  to  refine  it  within  the  boundary  layers.  We  further  assume  that 
derivatives  of  u  can  adequately  be  replaced  by  either  p(  t)  or  v(  a)  in  the  left 
or  right  boundary  layer,  respectively. 


This  problem  was  studied  by  Ascher  and  Weiss  (Ref.  2j  who  showed  that’ 
Eq.  (24)  could  be  approximately  satisfied  in  the  left  boundary  layer  by 
choosing  subinterval  lengths  as 

C  S( 1+| |u(ti_j) | | )  l/2k 

ti  *  ti-1  =  (--)(— -r-;— (26) 
o_  c) 1 w( t i_l ) j 1 


for  collocation  at  the  image  of  k  Gauss-Leger.dre  points  per  subinterval.  Here 
c  is  a  numerical  constant  and  a_  is  the  magnitude  of  the  largest  diagonal 
element  of  T_(X(0)).  A  similar  formula  can  be  obtained  for  selecting 
subinterval  lengths  in  the  right  boundary  layer. 

Starting  with  i  *  1  we  use  Eq.  (26)  to  generate  a  partition  until  we 
either  reach  t  =  1/2  or  a  point  where  a  subir.lerval  length  selected  by  Eq. 

(26)  is  larger  than  that  used  by  COLSYS  to  solve  the  reduced  problem.  We  then 
repeat  the  procedure  in  the  right  boundary  layer. 

We  have  written  a  computer  code  called  SPCOL  that  implements  the 
algorithms  that  are  described  in  this  section;  thus,  it  (i)  uses  COLSYS  to 
solve  the  reduced  problem,  (ii)  calculates  and  adds  appropriate  boundary  layer 
corrections  to  the  reduced  problem,  and  (iii)  (optionally)  suggests  a  mesh 
that  can  be  used  by  COLSYS  to  solve  the  complete  problem, 

4.  EXAMPLES.  In  order  to  appraise  the  performance  of  SPCOL  we  have 
applied  it  to  several  examples  involving  the  deformation  of  a  nonlinear 
elastic  beam  which  is  resting  on  a  nonlinear  elastic  foundation  and  is 
subjected  to  the  combined  action  of  a  horizontal  end  thrust  P  and  a  lateral 
load  p(x,t)  per  unit  length  (cf.  Figure  1).  This  problem  is  discussed  and 
analyzed  in  detail  in  Flaherty  and  O'Malley  (Ref.  6)  and  herein  we  only 
present  the  governing  equations,  which  in  diner.s ionless  form  are 


V 


X1 


•  • 

-  cos  x  ,  X2  -  Sin  X3  ,  X3  «  yi  (27a, b,c) 

eyj  -  -y2  ,  ev2  »  ( X^-p)  c°8  *3  -  Tyi  ,  (27d,e) 

where 

T  =  sec  X3  +  cv2  tan  X3  (27f) 

The  slow  variables  (xj,x2)  and  X3  represent  the  Cartesian  coordinates  and  the 
tangent  angle  of  a  material  particle  on  the  centerline  of  the  beam  that  was  at 
the  Cartesian  location  (t,0)  in  the  undeforraed  state.  The  fasL  variables  y\ 
and  y2  are  the  internal  bending  moment  and  transverse  shear  force, 
respectively  (cf.  Figure  1).  Finally,  the  snail  parameter  is 

e2  ■=  EI/PL2  ,  (28) 


where  El  is  the  flexural  rigidity  and  L  is  the  length  of  the  beam;  thus,  our 
beam  is  much  stronger  in  extension  than  it  is  in  bending. 

This  example  does  not  precisely  fit  out  hypotheses  since  the  axial  force 
T  is  a  function  of  the  fast  variable  y2  and,  thus,  G2'also  depends  on  y. 
However,  our  theory  and  methods  will  still  apply  as  long  as  y  remains  bounded 
and  |x3|  <  ir/2  as  e  *  0.  In  order  to  illustrate  the  diverse  behaviors  that 
can  occur  when  y  either  does  or  does  not  remain  bounded  as  e  +  0  we  present 
solutions  for  two  problems  both  having  X  =  p  =  1  and  which  differ  only  in 
their  boundary  conditions.  Some  additional  examples  are  presented  in  Flaherty 
and  O’Malley  (Refs.  6  and  8). 


In  our  first  example  we  take  the  boundary  conditions  as 


xi(0)  =  0  ,  -10x2(0)  +  y2(0)  =  0  ,  -X3(0)  +  10yi(0) 

10x2( 1)  +  y2(l)  -  0  ,  10x3(1)  +  y i ( 1 )  -  0 


0 

(29) 


These  supports  correspond  the  a  beam  that  is  almost  simply  supported  at  t  «  0 
and  almost  clamped  at  t  <*  1.  However,  perhaps  due  to  friction,  there  is  some 
coupling  between  lateral  and  rotational  effects  at  the  supports. 


As  we  shall  see,  y  remains  bounded  in  this  example  so  our  methods  are 
applicable.  The  orthogonal  matrix 


E(x(0) ) 


( l+a2)"l /2 


1  -a 
a  1 


(30a) 


a 


2 


where 


sec  X3(0) 


(30b) 


reduces 


G2(x( 0) ,0,0) 


(31) 


to  the  Schur  form  given  by  equation  (7)  ai  t  •  0  and  E^  will  reduce 
G2(x(1),1,0)  to  the  form  given  by  F.q.  (17)  at  t  »  1. 

We  solved  this  problem  in  two  ways:  (1)  using  COLSYS  to  solve  the 
complete  problem  (Eqs.  (27)  and  (29))  with  continuation  from  a  large  to  a 
small  value  of  e  and  (11)  using  our  code  S"COL  to  compute  an  initial 
asymptotic  approximation  and  to  recommend  a  nonuniform  mesh  and  using  this 
with  COLSYS  to  calculate  an  improved  solution.  All  calculations  were 
performed  in  double  precision  on  an  IBM  3033  coinput >r,  used  two  collocation 
points  per  subinterval,  and  set  the  error  tolerance  S  (cf.  Eq.  (24))  at  10-6 
for  slow  variables  and  10-3  for  fast  variables. 


Our  results  for  the  normalized  CP  times  and  the  number  of  subiniervals 
(NSUB)  that  are  either  used  by  COLSYS  or  recommended  by  SPCOL  are  shown  ir. 
Tables  1  and  2  for  continuation  in  e  and  our  methods,,  respectively.  Differ¬ 
ences  between  our  initial  asymptotic  approximation  and  the  final  solution 
obtained  by  COLSYS  are  shown  for  *3  and  yj  aL  t  *  0  and  1  in  Table  3.  We  see 
that  the  differences  decrease  like  0(e)  as  expected.  Dif'f erer;es  that  are 
recorded  as  zero  are  less  than  10-9.  Finally,  we  exhibit  solutions  for  x2 , 
x3 1  yi ,  and  y2  in  Figure  2. 

The  results  reported  in  Tables  1  and  2  need  some  additional  explanation. 
The  number  of  subiniervals  and  CP  times  used  with  continuation  depended  heav¬ 
ily  on  the  e  sequence  that  was  used.  The  results  in  Table  1  are  about  the 
best  insofar  as  they  gave  the  smallest  total  CP  time  for  the  sequence.  In 
addition,  COLSYS  relies  on  the  difference  between  solutions  that  are  computed 
on  two  different  partitions  in  order  to  estimate  local  errors.  Thus,  at  a 
minimum,  COLSYS  would  always  double  our  suggested  mesh.  This  is  apparent  in 
the  results  listed  under  the  heading  of  "COLSYS  Correction  No.  1"  in  Table  2. 
In  some  sense  these  results  are  encouraging  insofar  as  they  indicate  that  our 
mesh  selection  strategy  is  doing  about  as  well  as  it  can,  at  least  for 
e  <  10~2 .  However,  it  seems  that  fewer  points  should  be  necessary,  so  we 
tried  giving  COLSYS  an  initial  mesh  that  consisted  of  every  other  point  of  our 
recommended  mesh.  This  is  clearly  a  risky  strategy  since  collocation  at  the 
Gauss-Legendre  points  is  known  to  be  unstable  unless  the  mesh  is  sufficiently 
fine  in  the  boundary  layers  (cf.  Ascher  and  Weiss  (Ref.  2)).  Our  results 
using  this  are  reported  under  the  heading  of  "COLSYS  Correction  No.  2”  in 
Table  2.  Some  improvement  is  noted  for  e  >  1U-4;  however,  COLSYS  failed  to 
find  a  solution  (within  our  prescribed  limitations)  when  c  =  10-8. 

In  our  second  example  we  use  the  boundary  conditions 


xj(0)  =  0  ,  -x2(0)  +  ev2(0)  =  0  ,  -X3CO)  +  e2yi(0)  -  0 
x2(l)  +  ey2( 1 )  »  0  ,  X3(l)  +  EZyi(l)  -  0 


(32) 


TABLE  1.  NONLINEAR  ELASTICALLY  SUPPORTED  BEAM.  NUMBER  OF  SUBINTERVALS  (NSUB) 
AND  CP  TIMES  USED  TO  SOLVE  THE  PROBLEM  BY  COLSYS  WITH  CONTINUATION 
IN  E.  THE  TOTAL  CP  IS  THE  ACCUMULATED  TIME  FOR  THE  e  SEQUENCE. 


e 

NSUB 

CP 

Total  CP 

10"1 

80 

8.0 

8.0 

10-2 

78 

9.0 

17.0 

10"4 

78 

19.5 

36.5 

10"6 

156 

44.5 

81.0 

10"8 

100 

19.0 

100.0 

TABLE  2.  NONLINEAR  ELASTICALLY  SUPPORTED  BEAM.  NUMBER  OF  S'JBINTERVALS  (NSUB) 
AND  CP  TIMES  TO  SOLVE  THE  PROBLEM  BY  SPCOL  AND  OBTAIN  AN  IMPROVEMENT 
BY  COLSYS.  THE  CP  TIMES  FOR  SPCOL  INCLUDE  THE  TIME  TO  CALCULATE  THE 
REDUCED  SOLUTION  WHICH  WAS  4.8  TIME  UNITS.  CORRECTION  NO.  1  USES 
THE  MESH  THAT  WAS  RECOMMENDED  BY  SPCOL.  CORRECTION  NO.  2  USES  A 
MESH  THAT  IS  TWICE  AS  COARSE.  THE  TOTAL  CP  IS" THE  SUM  OF  THE  TIMES 
FOR  THE  SPCOL  AND  COLSYS  SOLUTIONS. 


SPCOL 

COLSYS 

Correction  No.  1 

COLSYS 

Correction  No.  2 

e 

Rec.  No. 
of  NSUB 

CP 

NSUB 

_ J 

CP 

Total  CP 

NSUB 

CP 

Total  CP 

10"1 

40 

bi 

'12.0 

lo  .9 

80 

12.1 

16.9 

10"2 

45 

12.0 

16.9 

78 

8.1 

12.9 

10"4 

54 

16.9 

21.8 

66 

9.2 

14.1 

10"8 

55 

H 

110 

17.5 

22.3 

Failed 

TABLE  3.  NONLINEAR  ELASTICALLY  SUPPORTED  BEAM.  DIFFERENCES  BETWEEN  SPCOL 
AND  COLSYS  SOLUTIONS,  WHERE  A(  )  :=  |(  )SpC0L  ~  (  )cOLSYsl 


£ 

Ax3(0) 

Ay2(0) 

Ax3(  1) 

Ay3(  1) 

10" 1 

3.3xl0"1 

5 . 1 x 1 0“ 2 

6.8x10" 1 

3.6xl0_1 

10"2 

2.8xl0~2 

6. 6x10" 3 

6 . 1  x  1 0"  2 

3.9xlO"2 

10"4 

2. 7xl0"4 

6.8x10" 5 

6.  lxlO"4 

3.9xl0"4 

10"8 

0 

1.3x10” 7 

0 

0 
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Figure  2.  Numerical  solution  of  elastically  supported  beam  with 
boundary  conditions  given  by  Equations  (29). 
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Figure  3.  Numerical  solution  of  elastically  supported  beam  with 

boundary  conditions  given  by  Equations  (32) .  Note  that 
Y’l  and  y0  are  multiplied  by  e. 
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